Cavity-water interface is polar 
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We present the results of numerical simulations of the electrostatics and dynamics of water hy- 
dration shells surrounding Kihara cavities given by a Lennard- Jones (LJ) layer at the surface of a 
hard-sphere cavity. The local dielectric response of the hydration layer substantially exceeds that 
of bulk water, with the magnitude of the dielectric constant peak in the shell increasing with the 
growing cavity size. The polar shell propagates into bulk water to approximately the cavity radius. 
The statistics of the electrostatic field produced by water inside the cavity follow linear response 
and approach the prediction of continuum electrostatics with increasing cavity size. 
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Nanoscale interfaces of polar liquids combine strong 
distortions of the liquid density profile with highly per- 
turbed long-range electrostatic correlations. The inter- 
facial density, and the related diffusional dynamics, are 
governed by short-range packing restrictions and are rel- 
atively short-ranged. Nevertheless, the density fluctua- 
tions of the interfacial region are critical for the long- 
range hydrophobic forces [l|, Q and the related weak 
dewetting of interfaces of non-polar solutes [1, 0] ■ In con- 
trast, electrostatic interactions, and the orientational cor- 
relations of multipolar moments, are long-ranged. They 
are assigned macroscopic length-scale in the Maxwell 
(continuum) electrostatics propagating the effect of par- 
tial charges at dielectric interfaces on the length-scale 
of the Coulomb potential. Whether this picture is cor- 
rect for polar liquids and how the surface polarization is 
screened by the mobile liquid dipoles remains an open 
question [5j , the resolution of which will define the limits 
of continuum electrostatics in application to nanoscale 
interfaces of molecular liquids. 

Water presents a particular challenge to the prob- 
lems of interfacial dynamics and thermodynamics since 
energetically strong hydrogen bonds add a short-range 
scale competing with long-range electrostatic forces. The 
properties of hydration layers surrounding nanoscale so- 
lutes indeed turn out to be unusual. Apart from ubiq- 
uitous hydrophobic interactions linked to the structure 
of the water interface [H, measurements of microscopic 
electrostatics of the protein/ water interface have shown 
some surprising results. Electrostatics on the microscopic 
scale is traditionally probed by the dynamic and static 
band-shifts of optical dyes [6] . The corresponding Stokes- 
shift dynamics at the protein/ water interface showed a 
long exponential decay absent for the free chromophores 
in solution. This observation has prompted the label of 
"biological water" for hydration layers of biopolymers Q ■ 
While the cause of this effect is still debated [8j,[9j, various 
extent of slowing of the collective Stokes shift dynamics 
has been universally observed at protein/ water interfaces 
Further, the hydration shells around proteins were 



found to carry high local polarity [lfjj, thus linking the 
slower dynamics to the structural reorganization of wa- 
ter in the form of a polarized cluster around proteins 

11]. Unfortunately, the problem of the protein hydra- 
tion is inseparable from the complex protein dynamics 

12]. Studies excluding this latter component are there- 
fore necessary to understand the electrostatics of hydra- 
tion layers when the solute size grows to the nanoscale. 
This is the goal of this report. 

Here we present extensive Molecular Dynamics (MD) 
simulations of the structure of water around non-polar 
solutes (cavities). In contrast to previous active research 
in this field dSEl, we ask here the following questions: 
(i) how polar is the interface? (ii) how far into the bulk 
does the polarity perturbation propagate? and (iii) how 
are the orientational dipolar dynamics of the hydration 
layers affected by the solute? The main result of this 
study is the observation of a significant increase of the 
local water polarity at the interface, with the region of 
enhanced polarity extending into the bulk to approxi- 
mately the cavity radius. 

The water nanoscale interface was modeled by insert- 
ing spherical solutes carrying a hard-sphere (HS) core 
surrounded by a Lennard- Jones (LJ) potential layer. The 
interactions with the SPC/E oxygen is then given by the 
Kihara solute-solvent potential 
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The LJ well has the width a = 3 A and the energy €lj 
for which two values were used: £lj = 0.65 kJ/mol equal 
to the LJ energy between oxygens of SPC/E water and 
£lj = 20 kJ/mol close to the energy of hydrogen bonds 
in bulk water. The cavity size was varied by changing the 
HS radius rns m the range 0-12 A. The number of waters 
in the simulation cell was varied to allow sufficiently large 
solvation layers, with 4053 and 11845 hydration waters 
used for the smallest and largest solutes, respectively. 
The trajectories for analysis were 5 ns long, following 
100-500 ps equilibration. Simulations were performed 
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FIG. 1. The second-order orientational order parameter "p\ of 
the water first shell vs the cavity radius R — ms + <?. The 
solid diamonds and triangles refer to cavities in water with 
£lj = 0.65 and 20 kj/mol, respectively. The open points 
refer to cavities in the fluid of dipolar hard spheres [f| with 
the reduced dipole moments (m*) 2 = /3m 2 /erf equal to 2.0 
(open squares) and 3.0 (open circles); m is the dipole moment 
and <T S is the hard-sphere diameter of the solvent. The inset 
shows the orientational parameter p{ of the first hydration 
layer. 



with cubic periodic boundary conditions, at 273 K and 
zero pressure, with Berendsen thermostat and barostat 
and a timestep of 2 fs. Ewald sums with tin foil boundary 
conditions were used for electrostatic interactions. 

A spherical solute induces a spherical symmetry 
breaking in an otherwise isotropic liquid. The 
orientational structure of the interface consistent 
with this imposed symmetry is characterized by the 
first- and second-order orientational order parame- 
ters: pi(r) = (N(r))- 1 (j2 rj <r*j ■ *%') and p 2 (r) = 

(2N(r))~ 1 (E ri <r[ 3 (^ ' ^jf ~ These parameters 
project the unit dipolar vectors rhj within the shell of 
radius r on the radial direction fj = rj/rj, N(r) is the 
number of waters within the shell. The first hydration 
layer is then defined as R < r < R + 1.5 A, R = rns + & 
and the corresponding order parameters are p[ and p\ . 

The dielectric constant of a cavity-water mixture can 
be defined in terms of the volume occupied by the cavity 
relative to the volume of water 14]. The solute-solvent 
response function describing the interfacial polarization 
can then be calculated by accounting for dipolar fluc- 
tuations accumulated within a radial water layer sur- 
rounding the cavity x(r) = /3((5M(r)) 2 )/(3V(r)); /3 = 
I/^bT) is the inverse temperature. The water dipole 
moment M(r) in this equation is taken over the radial 
shell between the spherical cavity and radius r extend- 
ing from the cavity center to the bulk, V(r) is the shell 
volume. The limit of infinite dilution yields the bulk di- 
electric susceptibility of water x = x(°°) = ( e — 1)/(4tt), 
where e is the water dielectric constant. 

The fluctuation susceptibility x( r ) determines the di- 
electric constant e(r) accumulated within the radial layer. 
For simulations employing periodic boundary conditions 
with tin-foil boundary around replicas of the simulation 
cell the connection between the simulated variance of the 



dipole moment and the dielectric constant is particularly 
simple 15]: e(r) = 1 + 47rx(r). The macroscopic dielec- 
tric constant of water is then e = e(oo). 

The orientational structure of polar liquids at inter- 
faces is strongly affected by short-range orientational cor- 
relations. Unsaturated hydrogen bonds of surface wa- 
ters produce preferential in-plane orientations of water's 
dipoles [16| as reflected by the first and second orien- 
tational order parameters (Fig. [1]). The first-order pa- 
rameter p{ is nearly zero pointing to no preferential ra- 
dial orientation (inset in Fig. [I}, while p\ is non-zero and 
negative, in accord with the preferential in-plane orienta- 
tion of the dipoles. This orientational pattern is specific 
for water and does not necessarily repeat itself in other 
polar liquid. For comparison, p\ °f hard-sphere dipoles 
at the surface of a spherical cavity [f| passes through a 
minimum (Fig. [T|). Orientational order first grows with 
increasing cavity size as frustrations of dipolar orienta- 
tions are released with the growing number of dipoles. 
However, further increase of the cavity size leads to a 
weak dewetting of the interface by the puling force of the 
liquid 0] with the resulting destruction of the interfacial 
orientational order. In contrast, water preserves its par- 
allel interfacial order, mostly determined by its hydrogen- 
bond network and not much effected by the strength of 
the solute-solvent LJ attraction (Fig. [1]). 

The function e(r) calculated for shells around cavities 
is compared to the same function calculated from shells 
around Lorentz's virtual cavity [13] (water molecules be- 
tween radii R and r) taken from configurations of pure 
water without cavity inserted. The difference of the two 
functions shows a sharp peak (Fig. [2h) pointing to an ef- 
fectively higher polarity of hydration shells around cavi- 
ties compared to shells in bulk water. 

In order to distinguish between the orientational and 
density origins of the peak in the dielectric constant, 
we have plotted in Figs. EJb,c) xi r ) defined as above 
in comparison with the susceptibility normalized to 
the number of waters in the shell N(r): xn(t) ~ 
(3p{(SM.(r)) 2 ) /(3N(r)), where p is the number density of 
bulk water. This comparison shows that the origin of Ae 
is a composite effect of changes in both the local density 
and orientational structure. Even though a significant 
part of Ae comes from the increased density in the first 
solvation layer, particularly at the large solute-solvent LJ 
attraction (Fig. |2](c)), the effect cannot be cast in terms 
of N{r) only. 

Given the long range of the interfacial orientational 
order one wonders how the dynamics of hydration layers 
are affected. We have looked at several correlation func- 
tions. x'(t) = PdSM'it) ■ (5M / (0))/(3U / ) is the time 
self-correlation function of the dipole moment of the first 
solvation layer and x(r, t) is a similar correlation function 
extended to a layer within the radius r from the cavity's 
center. In addition, we have calculated the correlation 
function Cs(t) = (<5E s (r, t) -<5E s (r, 0)) of the electric field 
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FIG. 2. Dielectric constant of the hydration layer relative to 
the dielectric constant of the same layer around a virtual cav- 
ity for three cavity sizes indicated in the plot (a). The solid 
and dashed lines refer to elj = 20 and 0.65 kj/mol, respec- 
tively. The response functions defined through the volume of 
the shell ("V", solid lines) and through the number of shell 
waters ("N", dashed lines) are shown in (b) for elj = 0.65 
kJ/mol and in (c) for e L j = 20 kj/mol; R = 7.5 A. xM for 
the virtual cavity is marked as "virt." 
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FIG. 3. Exponential relaxation time of x 1 if) (closed dia- 
monds), Cs(r, t) (open circles), and x( r i^) (open squares); 
£lj = 20 kJ/mol. Also shown is the exponential relaxation 
time of the self-correlation function of the unit vector e 1 (t) 
representing the dipole moment of the first-shell waters. Dif- 
ferent cavity sizes for the r-dependent relaxation times are 
indicated in the plot. The filled circles refer to R = 12 A and 
6l.t = 0.65 kj/mol and the horizontal dotted line indicates the 
Debye relaxation time td of pure SPC/E water. The dashed 
lines in the plot are connecting the points. 



E s (r, /) produced at the cavity's center by the waters 
within the r-shell. This latter correlation function repre- 
sents the Stokes shift dynamics of dipolar chromophores 
@, HI ■ All correlation functions were fitted to a sum of 
a ballistic Gaussian decay and an exponential tail @. 
Figure [3J presents the compilation of the results for the 
exponential relaxation time te ■ 

The main observation from the dynamics calculations 
is a significant growth of the exponential relaxation time 
of the first solvation layer with the cavity size (filled dia- 



monds in Fig- |3j) . Consistent with the results for Ae(r), 
this dynamics perturbation propagates into the bulk to 
at least the distance of the cavity radius. The expo- 
nential decay time of CE(r, /) retraces the correspond- 
ing time from x{ r i t)- The dipolar dynamics of the first 
solvation layer are dominated by rotations of the dipole 
moment M 7 , instead of its magnitude fluctuations, as is 
seen from the self-correlation function of the unit vector 
e T (t) = M J (i)/M 7 (i) (crosses in Fig. [3]). This observa- 
tion points to a high level of orientational cooperativity 
in the first solvation layer, which does not decorrelate by 
individual dipole rotations and instead rotates slowly as 
a correlated dipolar domain. This dynamical slowing is 
however seen only for the larger LJ attraction, 6lj = 20 
kJ/mol, and no effect of the solute on the dynamics is 
observed when the solute-solvent LJ potential is similar 
to that in water, £lj = 0.65 kJ/mol (closed points in Fig. 
|3j . This result might help to explain the conflicting liter- 
ature [ItJ on the subject of the surface-induced alteration 
of the liquid dynamics. The outcome seems to be con- 
trolled by the strength of the solute-solvent interactions 
and thus the surface composition. 

Our results partially support Onsager's concept of "in- 
verted snowball" dynamics (l8| . It stipulates that solva- 
tion dynamics are slower close to a newly created charge 
compared to more distant layers, ranging between the 
one-particle (slow) orientational diffusion and the dielec- 
tric (fast) relaxation of the bulk. The data in Fig. [3] in- 
deed show a speedup of dipolar relaxation into the bulk. 
However, this effect strongly depends on both the cavity 
size and the solute-solvent LJ interaction. It is expected 
to be essentially absent for typical optical probes [6| con- 
sistent in size with the smallest cavity studied here. Fur- 
ther, even for larger cavities, the slow dynamics of the 
closest solvation layers are almost lost when the hydra- 
tion layer is grown to the boundaries of the simulation 
cell. At that point, the relaxation time becomes the De- 
bye relaxation time of bulk water (Fig. [3J. This obser- 
vation implies that dipolar optical probes placed inside 
the cavity 7-9] will not pick up the slowing of the closest 
hydration shells and instead will average the effect out by 
the electric field contributions from more distant layers. 

The statistics of dipolar interfacial fluctuation can be 
probed by the chemical potential of electrostatic solva- 
tion, i.e. the free energy of interaction of the charges in- 
side the cavity with the surrounding water solvent. We 
found that the dipolar field M(r) is Gaussian and thus 
the linear response approximation should be applicable. 
The solvation chemical potential of a charge [i q or a 
dipole Hd can then be found from the variance of the elec- 
trostatic potential 4> s or the electric field E s produced by 
the solvent at the position of the corresponding multipolc 
0: M? = -/3(g ) 2 <W> s ) 2 )A fi d = -/3(m ) 2 ((5E s ) 2 )/6. 
Here, go and tuq are, correspondingly, the charge and 
point dipole within the cavity. The averages in these re- 
lations do not depend, in linear response, on whether the 
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FIG. 4. Chemical potential fi" d = fidrn^/R 3 of solvating 
the point dipole mo at the center of the cavity of radius R. 
The solid points = 20 kj/mol, circles and 0.65 kj/mol, 
squares) and obtained from the variance of the water electric 
field inside the empty cavity, fid oc ((SE S ) 2 }. The open trian- 
gles are from the thermodynamic integration of the average 
electric field (E s ) produced by dipoles of increasing magni- 
tude (0 < mo < 10 D) positioned at the cavity center. The 
average field (eV/D) is a linear function of the dipole magni- 
tude mo (inset, R = 12 A, elj = 0.65 kj/mol). The dashed 
line in the inset is the prediction based on the field variance 
inside the empty cavity (not the best fit). The dotted hori- 
zontal line is the result of continuum electrostatics given by 
the Onsager relation ("O"). 



corresponding multipoles are actually present inside the 
cavity [l9j. They can therefore be calculated from our 
simulations with empty cavities providing insights into 
how these solvation free energies scale with the size of 
the solute and whether the limit of continuum electro- 
statics is reached. 

We found noticeable effects of the size of the sim- 
ulation cell on {(5(f> s ) 2 ) and much smaller size effects 
on ((<5E S ) 2 ). We have therefore chosen to look at the 
statistics of the field fluctuations and the corresponding 
chemical potential of dipole solvation. These results are 
summarized in Fig. [4] showing fi* d = /i^i? 3 / '( TO o) 2 • This 
dimensionless parameter is expected to approach, with 
growing cavity size, the size-independent limit of contin- 
uum electrostatics, given by the Onsager equation [l4[ 
fid = (e — l)/(2e + 1) — 0.5. The values of fid, al- 
though noticeably higher at intermediate sizes, indeed 
seem to approach this limit. The open points in Fig. 0] 
are obtained by thermodynamic integration of average 
energies of point dipoles placed at the cavity center. A 
good agreement between fid from the field variance and 
from the thermodynamic integration, as well as the lin- 
ear dependence of (E s ) vs too (inset in Fig. 2]), testifies 
to the validity of the linear response approximation. 

The chemical potential fid is not strongly affected by 
the strength of the solute-solvent LJ attraction. The 
range of clj values studied here probably covers most 
of situations of practical interest. Deviations of the elec- 
tric field variance from the area between the two curves 



shown in Fig. [4] might therefore be used to identify nonlin- 
ear solvation typically associated with hydrogen bonding 
between water and the solute [2(| . 

The results obtained here must have significant impli- 
cations for self-assembly of nano-sized objects and bio- 
logical activity of hydrated biopolymers. Polar solvation 
layers around solutes are expected to screen the inside 
charges. In the crowded environment of a living cell 21] 
this screening will reduce interactions between multipo- 
lar solutes. The high-polarity layer is also characterized 
by slower dipolar solvation. However, this effect is is not 
picked up by the dynamics of dipolar probes placed in- 
side the cavity. A slow relaxation component observed 
in optical time-resolved spectra [7| therefore needs to be 
assigned to protein motions pushing the hydration layers 
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